Systematic characterization of thermodynamic and dynamical phase behavior in 

systems with short-ranged attraction 
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In this paper we demonstrate the feasibihty and utihty of an augmented version of the Gibbs 
ensemble Monte Carlo method for computing the phase behavior of systems with strong, extremely 
short-ranged attractions. For generic potential shapes, this approach allows for the investigation 
of narrower attractive widths than those previously reported. Direct comparison to previous self- 
consistent Ornstein-Zernike approximation calculations are made. A preliminary investigation of 
out-of-equilibrium behavior is also performed. Our results suggest that the recent observations of 
stable cluster phases in systems without long-ranged repulsions are intimately related to gas-crystal 
and metastable gas-liquid phase separation. 
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I. INTRODUCTION 

Colloidal systems with short-ranged interactions serve 
as a minimal model of complex fluids used in a variety 
of technological applications [1]. While such systems are 
important in the applied realm, their intrinsic phase be- 
havior is of great fundamental interest. Upon adding 
polymer to an otherwise uniform colloidal^su^spension, an 
entropic depletion interaction is induced [2, 3]. The range 
and strength of this attractive interaction may be con- 
trolled by the length and concentration of the added poly- 
mer. Thus, exquisite control may be experimentally ex- 
ercised over such systems, which allows for an exhaustive 
exploration of the temperature-density phase diagram for 
nearly any range of attractive interaction 0, S Q • 

At high volume fractions, depletion attractions have 
been shown to induce an inverse melting behavior, such 
that the viscosity of a hard-sphere suspension close to its 
colloidal glass transition is significantly reduced [Tl, [8|. 
This suspension may be revitrified by addition of still 
more polymer. Thus, two glassy phases appear to exist 
at the same volume fraction: one induced by the repul- 
sions and one induced by strong attractions Q. It has 
been speculated that the attractive glass state at high 
volume fractions can be continuously connected to a gel 
state at lower volume fractions [9]. A major difficulty 
with the connection between attractive glass and colloidal 
gel is that phase separation often intervenes [10, nu\ - 
Phase separation into colloid rich and colloid poor re- 
gions may become anomalously slow if the density of 
the colloid rich region is close to the density of the uni- 
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form attractive glass at high volume fractions. This pro- 
cess leads to weak colloidal gels whose connection with 
near-equilibrium high volume fraction attractive glasses 
is complicated by the interplay of vitrification and phase 
separation [l2|,[l3- Thus, the added dimension of out-of- 
equilibrium behavior may greatly increase the complex- 
ity of the various phases that may be observed in such 
systems. 

In addition to colloidal gels induced by phase separa- 
tion, various cluster phases have been recently experi- 
mentally observed [3, |l5|, |l^ . If charge resides on the 
colloids, then equilibrium microphase clusters may form 
due to the competition between short-ranged attraction 
and long-ranged charge repulsion [17, 18, 19]. Interest- 
ingly, recent experiments that utilize sufficient salt to 
screen the charge on the colloids still show the existence 
of relatively stable large clusters. The existence of such 
clusters even in the absence of repulsion has been sug- 
gested theoretically [2^. Lu et al [15] have observed 
large, relatively compact clusters at low volume fraction 
and low (~ 1 — 2kBT) attraction strength. Sedgwick et 
al. [16] have observed a "bead phase" of clusters that 
appears to exist only along a portion of the metastable 
gas-liquid binodal. 

Clearly, the investigation of the nature of gel and clus- 
ter phases via computer simulation is complicated by 
the need for a precise characterization of the equilib- 
rium phase diagram. Qualitatively, the role of short- 
ranged attractions in widening the gas-solid coexistence 
gap is well known [2l|, [22| • The characterization of vari- 
ous phase boundaries for general short-ranged potentials, 
however, has been performed by integral equation meth- 
ods [2311 by a mix of Monte Carlo and analytic expan- 
sions [23, z5], on systems with a moderately short at- 
traction range [26, 27, 28, 29, 30, 31, 32, 33, 34, 35, ^, 
or short-ranged attractive systems with somewhat arti- 
ficial forms that do not lend themselves to dynamical 



studies [371, |38] . In this work, we will combine several ex- 
isting computational methodologies to produce a direct 
Monte Carlo approach that is capable of yielding essen- 
tially exact phase behavior over a wide range of the ther- 
modynamic parameter space, for systems with extremely 
short-ranged attractive interactions. 

In particular, we demonstrate the general applica- 
bility of a form of the Gibbs ensemble Monte Carlo 
(GEMC) [39| approach augmented with improved sam- 
pling techniques for the calculation of phase behavior 
in systems that mimic colloidal suspensions with short- 
ranged attractions. In addition, our approach allows for 
the direct sampling of configurations without additional 
a priori knowledge when the phase behavior is com- 
plex, such as in systems that exhibit microphase sepa- 
ration [^. We use the results gleaned from this imple- 
mentation of GEMC to compare with several techniques, 
including the approximate but powerful self-consistent 
Ornstein-Zernike approximation (SCOZA) approach [23] 
and exact approaches such as thermodynamic integra- 
tion p[]|. With these results in hand, we make a pre- 
liminary study of the nonequilibrium phase behavior of 
systems presumably similar to those studied by Lu et 
al. [15| and Sedgwick et al. [16] in the regimes where 
cluster phases might be expected. Our paper is orga- 
nized as follows. In Sec. [Ill we discuss the application of 
our GEMC methodology to systems with short-ranged 
attractions, and we outline the systems that are studied 
in this work. In Sec. Illli we compare the results of our 
GEMC to previously published results. In Sec. [IVl we 
study various aspects of the out-of-equilibrium behavior 
of systems similar to those published in recent experi- 
ments. Our studies are facilitated by the precise knowl- 
edge of the phase diagram afforded by the Monte Carlo 
method developed in this work. In Sec.[Vlwe conclude. 



II. SIMULATION TECHNIQUES 

Various potential shapes have been suggested to cap- 
ture the phenomenology of short-ranged attractive sys- 
tems. Two commonly used forms are the generalized 
Lennard- Jones (LJ) potential 
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with n > 12 [42 
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and the hard core plus attractive Yukawa 
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with b > 6.05 [23|. Both are used in the current work. 
Results are reported in reduced units. Temperatures are 
in units of the well depth e, distances are in units of the 
particle diameter <j, and time is rescaled by {e/ma'^y^'^ 
where m is the mass of the particles. 
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FIG. 1: GEMC 6 = 30 gas density for the gas-crystal coex- 
istence at T — 0.33 for different crystal slab dimensions N. 
Within error bars, the bulk limit (within the GEMC error 
interval) is attained when the crystal has > 500 particles. 



GEMC with single particle displacement and exchange 
as well as volume exchange is used for the metastable 
gas-liquid equilibrium for two boxes of 256 particles (500 
for 6 = 60 and LJ potentials). Finite size effects were 
checked by comparing the results to that of a square- 
well potential with a width of 0.25cr [23]. As the work 
was being completed, a paper with similar methodology 
(applied to square well fluids) appeared, and thus more 
details can be found there [29|. In the current work, 50% 
of the moves are particle swaps, 49.5% are particle dis- 
placements, and the rest are volume interchanges, while 
equilibrium and production runs involve at least 10^, and 
sometimes reaching 10^ Monte Carlo cycles. When solid 
nucleated in the liquid box, only the intermediate results 
are retained for analysis. 

For the gas-crystal equilibrium, the methodology de- 
veloped by Chen and Siepmann is followed [43, H] . This 
greatly increases the efficiency of GEMC in this regime. 
The boxes are configured such that one box contains only 
vapor and the other a slab of solid surrounded by vapor. 
In addition to basic GEMC moves (47% particle displace- 
ments, 40% particle swaps, 3% symmetric and asymmet- 
ric volume exchanges) aggregation-volume bias [43] and 
its generalization to exchanges between two boxes (5% 
each) are performed. For these simulations initial sys- 
tem sizes are 256 particles for the gas phase, 864 for the 
solid slab, and again a minimum of 10^ Monte Carlo cy- 
cles are performed for equilibration and production. 

All dynamical results are performed with the Lennard- 
Jones potential using standard molecular dynamics inte- 
gration for systems of size N = 5324 and with an in- 
tegration step At = 0.001. Cooling is done by velocity 
resampling every 10^ steps starting from a liquid configu- 
ration. We do not expect a strong influence of the type of 
dynamics on the qualitative results [45|, so further stud- 
ies with other dynamical protocols are not performed. 
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FIG. 2: (Color online) Comparison of GEMC results (green squares) with SCOZA [23] (solid black lines) phase diagram 
calculations for the hard sphere Yukawa potential with b = 100, 60, 30, and 7.5 for panel (a), (b), (c), and (d), respectively. 
The gas-crystal binodal for a given density is found at higher temperatures than the gas-Uquid binodal in all systems studied. 
The dashed line indicates the crystal close-packing limit. For b = 100, results from Ref. [41] (dash-dotted line) are also included. 
For b = 100 and 60, only the GEMC gas-crystal binodal is shown. The black dots are the critical point predicted by extended 
corresponding state, which are Tc = 0.1711, 0.1725, 0.231, and 0.382, respectively. The line between the squares is a guide for 
the eye in the gas-crystal case; refer to the text for the gas-Uquid case. The error bars due to density fluctuations are smaller 
than the symbol sizes. 



Finite-size effects 



An obvious source or error exists in the approach we 
have employed to determine the gas-crystal coexistence 
line. In particular, the standard use of GEMC for gas- 
crystal coexistence would utilize a box of pure solid and 
a box of pure vapor with particle exchanges between the 
two boxes. In such a configuration, exchanges between 
boxes would be prohibitively infrequent. We have in- 
stead used the Chen and Siepmann approach where the 
solid box is replaced by a box containing a solid slab sur- 
rounded by vapor [4J]. From the microscopic perspec- 
tive, several possible finite-size effects arise. First, the 
crystal block might show anisotropy in its vapor pressure 
depending on the selection of the grain surface exposed 
to the gas. However, this is expected to be very small 
for spherically symmetric potential shapes. Second, at 
high densities, crystal surface fluctuations or other im- 
age effects may become significant. This was avoided by 
working at relatively low gas densities and by using a 
box of sufficient dimensions to eliminate periodic arti- 
facts. From the thermodynamic standpoint, the Chen- 



Siepmann procedure can lead to biased results due to 
the fact that the free energy of the solid box now has 
a contribution from the surface. In principle, this error 
should scale as ~ 2(Js/L where as is the surface tension 
and L the crystal size. In previous studies, Chen and 
Siepmann have shown that their results are in quanti- 
tative agreement with the known gas-crystal coexistence 
behavior of several simple systems [4^ . While they did 
not perform a study of finite size effects, their results sug- 
gest that, in at least some instances, it is indeed possible 
to obtain quantitative results efficiently with a modified 
GEMC approach that allows for metastable equilibrium 
via "coexistence" between a gas on the one hand and a 
gas-crystal mixture on the other. However, the conver- 
gence as a function of system size should be performed at 
each temperature and for each system, in order to ensure 
accurate results within the modified GEMC approach of 
Chen and Siepmann. Here, we carry out such a study. 

In Fig. [1] we show the gas density obtained for b = 30 
and T = 0.33 for a variety of different system sizes. Af- 
ter some noticeable damped oscillatory behavior for small 
system sizes, the value remains within the GEMC error 



width. We consider the density to be converged when 
the result has saturated within a band of densities whose 
width is that of the error bars of the GEMC procedure 
itself for a given size of the solid slab. The nontrivial con- 
vergence behavior hints at a possible lattice periodicity 
effect (e.g., the period of the behavior might be related 
to the number of particles in a layer of solid) , but a more 
detailed analysis is left for a later study. The results pre- 
sented in Fig. [1] are representative of the results found for 
other values of b and T that are presented in this work. 
Thus, we are confident of the accuracy of the results 
for the gas-crystal lines presented in Figs. [2] and [3l While 
it is true that the method used here for gas-crystal coex- 
istence can induce systematic error, it appears that such 
errors can be kept to a controllable size without undue 
computational effort, at least in the systems we have in- 
vestigated. It should also be noted that while approaches 
such as thermodynamic integration do not suffer from 
the potential systematic bias discussed in this subsection, 
they still may lead to inaccurate results if a good refer- 
ence system is not used (as could the case for extremely 
short-ranged attractive potentials) or if the systems uti- 
lized are too small. This will be discussed further in the 
next section. 



III. COMPARISON WITH SCOZA 




FIG. 3: (Color online) GEMC phase diagram (green squares) 
for the LJ potential with n = 50 is shown superimposed with 
shifted SCOZA (solid black lines) results for b — 30. The lines 
between the squares are a guide for the eye. The SCOZA 
spinodal line is indicated by a dash-dotted line. See text for 
details on the fit. The error bars from density fluctuations 
are smaller than the symbol sizes. 

Various aspects of the thermodynamic properties of 
systems with short-ranged attractions have been studied 
in the last decade [22|,M 0, S, S, [s^, [Ml, S. One 
of the interesting features of such systems is the possi- 
bility of metastable gas-liquid coexistence below the gas- 
crystal equilibrium region. However, no verification of 
integral equation predictions, other than in the Baxter 





FIG. 4: (Color online) The system at point B' in Fig. [3] once 
nucleated at (a) early times and (b) later times. The details 
of the artificial seeding process are discussed in the text. 



limit of an infinitely narrow potential, has been made in 
extremely short-ranged attractive systems [38]. SCOZA 
is generally viewed as being among the most accurate in- 
tegral equation methods available [46, 47]. Predictions 
of phase diagrams for short-ranged potentials with hard- 
core Yukawa have been made using SCOZA and SCOZA- 
based hybrid methods [23]. However, these results have 
only been compared with simulation results slightly over 
the limit of metastability with b = 7 [22|, [23, [3l| for the 
gas-liquid binodal and up to 6 = 9 for the gas- solid bin- 
odal [22]. The Baxter limit of infinitely small interaction 
range has also been reported via simulation [38] , but only 
for the gas-liquid binodal. For colloidal depletion inter- 
actions, ranges of the order of 0.25(T to O.Olcr are of more 
interest, corresponding to 6 = 20 — 100 for the hard-core 
attractive Yukawa. Simulated phase diagrams for sys- 
tems with an interaction range as low as 0.1 5cr for square 
well [25, 29], and > 0.08cr for Derjaguin-Landau-Verwey- 
Overbeek (DLVO) [33|, LJ-based [26, 30, 32], Asakura- 
Oosawa pair potential [SJ, [35|, [3^ have been reported, 
but none of these can easily be quantitatively compared 
with the SCOZA predictions, other than the cases pre- 
sented in Ref . [41| . 

Calculated phase diagrams for the hard core attractive 
Yukawa potential are shown in Fig. [2l The metastable 
gas-liquid density-temperature binodal was fitted using 
\p - Pel = A\T - Tc\ - B\T - Tcl^, where p = 0.3258 [48.] 
and the critical temperature Tc used in the fitting was 
extracted from the result of the Baxter limit [37] and 
mapped to this system by equating second virial coeffi- 
cients as done for similar systems [45] . Due to the flatness 
of the curve, critical densities pc are harder to pinpoint. 
A value of pc ~ 0.50(5) appears to fit all of our simu- 
lated systems. No unique set of A and B parameters fit 
the binodal through the entire temperature range of our 



simulations, but a rather reasonable agreement was ob- 
tained nonetheless for fixed A and B. The results for the 
hard-core Yukawa potential for 6 = 7 and 6 = 25 for Tc 
and pc are consistent with those of Ref. [41] for nearby 
values of b. 

For the gas-solid binodal, a simple spline fit was made 
to the gas branch. Densities on the solid side are indis- 
tinguishable from that of close packing within measured 
precision. Much larger system sizes would be required 
to obtain a sufficiently precise measurement. This leads 
us to suspect the solid phase density results obtained by 
a recent paper are due to an erroneous density calcula- 
tion scheme for these systems [29]. However, this is not 
the part of the phase diagram that is of interest here, 
and we will thus satisfy ourselves with a straight line at 
the close-packing density. Finally, the range of densities 
available for simulations did not allow for the study of 
solid- solid coexistence, known to take place in systems 
with an extremely short range of interaction [49|, \bO] . 

Comparing the phase diagrams of two different poten- 
tial shapes with a similar range of attraction (~ O.la), 
namely the LJ with n = 50 and the hard core attractive 
Yukawa potential with 6 = 30, we observe that the inter- 
action range is a good measure of gas-liquid phase separa- 
tion metastability: the gap between the metastable gas- 
liquid critical point and the gas-solid binodal is T ?^ 0.15 
in both cases. After shifting the phase diagram so as 
to match critical points, SCOZA predictions are seen to 
match the LJ simulation results just as well in Fig. [3l 
This is in the spirit of the extended corresponding state 
principle discussed above. 

Comparison of GEMC simulation with SCOZA shows 
rather good agreement for the gas-solid binodal, espe- 
cially for the wider interaction ranges of b — 7.5 and 
6 = 30. For the narrower ranges of 6 = 60 and 100, we 
reach densities beyond the range numerically attainable 
for SCOZA calculations [23]. Also, though the curves 
have similar shapes, the simulation results are shifted to 
lower temperatures with respect to the SCOZA predic- 
tions. On the high density side, phase boundaries would 
only be accessible through the GEMC methodology using 
much larger systems due to the fact that small blocks of 
crystal tend to melt completely at high vapor pressures. 
The pathological SCOZA predictions of solids beyond the 
close-packing density indicates the limitation of such an 
approach in high density regimes. 

GEMC metastable gas-liquid binodals for 6 = 7.5 and 
30 also agree well with SCOZA predictions, especially for 
the low density branch. It is difficult to reach lower tem- 
peratures as solid nucleation becomes facile in metastable 
finite-sized systems. This limits the range of temper- 
atures where the metastable binodal data can be ob- 
tained. Furthermore, this approach breaks down for nar- 
rower interaction ranges when the ffatness of the bin- 
odal does not allow for a sufficient separation of Monte 
Carlo time scales before the crystal nucleates. Obtaining 
spinodal curves by simulation would require precise dy- 
namical measurement upon cooling and as its meaning 



is loosely defined in any case [23], no direct calculation 
of the spinodal was attempted. For further reference the 
SCOZA spinodal was included in the LJ phase diagram 
in Fig. [3] [23]. 

Finally, we conclude this section with a brief discus- 
sion of the comparison of our results for a very narrow 
attraction range {b = 100) with that of thermodynamic 
integration techniques [41|. In principle, thermodynamic 
integration is the preferred method, since it does not have 
a systematic bias such as, for example, that discussed in 
Sec. Ill A[ However, in some important cases the meth- 
ods presented here have some advantages over this direct 
approach. In particular thermodynamic integration may 
be difficult to perform if the repulsive portion of the po- 
tential is not of the hardcore variety or if a suitable refer- 
ence system cannot be found. More importantly, in cases 
where the morphological characteristics of the phases are 
not known (such as in complex microphases or domain- 
forming systems [40]) it is of great advantage to have a 
direct Monte Carlo approach that does not need to make 
use of presupposed information. In Fig. [2^ we compare 
the thermodynamic integration results of Ref. [41] to that 
of our GEMC approach. While the results of Ref. [41] are 
consistent with our results for smaller values of 6, here, 
a difference of about 10% can be seen in the gas-crystal 
coexistence temperatures at low densities. We believe, 
based on an analysis of the type given in Sec. Ill Al that 
our results should be quantitatively accurate. Interest- 
ingly, a similar discrepancy of the same magnitude and 
direction exists between the location of the gas-liquid bin- 
odal line and the location of the critical point as found by 
the extended corresponding state predictions. It is pos- 
sible that the small system size {N = 108) used in the 
thermodynamic integration study is one source of these 
discrepancies. It would be of some interest to examine 
these issues further in a future study. 



IV. NONEQUILIBRIUM KINETICS IN 
VARIOUS REGIONS OF THE PHASE DIAGRAM 



The precise location of various phase boundaries has 
important implications for the dynamics that might be 
observed upon quenching a homogeneous system into dif- 
ferent temperature and density regions of the phase di- 
agram. Interesting out-of-equilibrium phenomena have 
been observed experimentally depending on the depth of 
the quench, the initial density, and the quench rate. Us- 
ing precise characterization of the phase diagram for the 
LJ (n = 50) system of Fig. [3l we use direct molecular 
dynamics simulation to investigate some aspect of this 
kinetic behavior. For clarity, we shall refer to the vari- 
ous labeled points in Fig. [3l Points A-E are located at a 
volume fraction (j) = np/G = 15% and at T = 0.50, 0.35, 
0.275, 0.20, and 0.10 from A to E, respectively. Points B' 
and C refer to volume fraction (j) = 6.9% and (p = 4.7% 
and temperature T = 0.3 and T = 0.26, respectively. 
The behavior of the system at point A, as expected, is 






(a) 



(b) 



(c) 



FIG. 5: (Color online) Progression of phase separation at T = 0.275 (point C in Fig. [3|) from homogeneous liquid to eventual 
gas-crystal phase coexistence via the spinodal decomposition due to crossing the metastable gas-liquid binodal. The system is 
at = 15% and is here shown at t = 0, 300, and 2000. For clarity monomers are represented by smaller spheres. 



that of a homogeneous gas as verified by direct dynam- 
ical simulation. No further mention of the behavior in 
this region of the phase diagram will be made. 



A. Gas-crystal phase coexistence 

Cluster phases in colloidal systems may exist for a va- 
riety of reasons [51, 52, 53, 54]. Perhaps the most ubiqui- 
tous reason is the presence of both short-ranged attrac- 
tions and long-ranged repulsion. This repulsion is most 
commonly thought to arise from excess charge residing 
on the colloidal particles. Recent experiments, however, 
have suggested that even when charge repulsion is highly 
screened, compact clusters may appear and evolve in a 
rather stable fashion. Lu et al. have recently studied, 
via confocal microscopy, the evolution of attractive col- 
loidal suspensions in a high salt environment [15]. They 
find evidence for stable, compact clusters of more than 
500 particles. This behavior is most prominent at rather 
weak attraction strength {U ~ 2k bT) and an interme- 
diate range of attraction (~ 0.15cr). Sedgwick et a/., 
in screened lysozyme solution, find evidence for a com- 
pact "bead" phase in a different region of the phase di- 
agram [iQ]. The bead phase of Sedgwick et al. appears 
to exist only close to the metastable gas-liquid binodal. 
Therefore, the ability to locate precise boundaries in the 
equilibrium phase diagram is a fundamental prerequisite 
for the study of such kinetic behavior. Below, we un- 
dertake a preliminary study of nonequilibrium behavior 
in the regions of the phase diagram where such cluster 
phases have been found. 

In the study of Lu et a/., the attraction strength where 
a cluster phase is observed would appear to be in the 
broad gas-crystal coexistence region [15]. Further sup- 
port for this may be found in the degree of crvstallinity 
seen in the clusters imaged by Lu et al. [GO]- In our 
simulations, we find that, by inserting a small face- 
centered-cubic nucleus, the subsequent crystal growth is 



rather facile and reaches equilibrium with the gas phase 
rapidly. The behavior occurs essentially throughout the 
gas-crystal coexistence region above the metastable gas- 
liquid binodal. The idea that nucleation and growth 
of clusters can be self-limiting, leading to a metastable 
cluster fluid phase, was theoretically put forth by Kroy, 
Cates, and Poon [2^. Lu et al. propose that the rear- 
rangement of particles at the surface and on the interior 
of a cluster may occur on a different time scale than the 
time scale on which clusters diffuse away from each other, 
rendering clusters long-lived. This notion is similar to 
that of Ref. 0. 

To test this kind of hypothesis we study larger system 
than that possible by GEMC. We artificially nucleate the 
solid by carving out irregular portions of the equilibrium 
crystal (~ 2000 particles each) allowing the surface of 
these clusters to come to equilibrium with the vapor un- 
der the appropriate thermodynamic conditions, and then 
allowing the clusters evolve via molecular dynamics. We 
have implemented this procedure at points B and B' of 
the phase diagram in Fig. [3l While this procedure is a 
crude mimic of the kinetics of gas-crystal nucleation, it 
does provide a useful test of self-limited cluster growth 
as a stabilizing mechanism of the cluster phase observed 
by Lu et al. [15|. We find that smaller compact clus- 
ters merge without much impediment upon collision with 
each other. This process is depicted in Fig. [H No tra- 
jectories showed cluster dissociation. It should be noted 
that the same behavior is observed at higher tempera- 
tures as well, where the initial cluster surface morpholo- 
gies are rougher due to a higher equilibrium vapor pres- 
sure. We did not find any evidence of a strong cluster size 
dependence for this process. Thus, while we did not make 
a systematic study of the effects of quench rate, system 
size, or polydispersity, one possibility for the cluster fluid 
phase observed by Lu et al. is simply a slowly evolving 
solid nucleation and growth process [60]. Interestingly, 
Sedgwick et al. find that the crystal phase is observed 
for a wide range of volume fractions and temperatures 



effectively in the gap between gas-solid and metastable 
gas-liquid coexistence. This is consistent with our find- 
ings in the analogous region of the phase diagram. 

A final possibility that should be mentioned with re- 
gard to the cluster phase observed by Lu et al. is accumu- 
lation repulsion [55|, l56|. Even in the absence of charge, 
it is possible that the neutral added polymer may induce 
an effective repulsion between colloidal spheres, due to 
its enhanced concentration at the surface of the colloids. 
It would appear unlikely that this repulsion is sufficient 
to render clusters thermodynamically stable in the sys- 
tem studied by Lu et al. However, it is possible that 
this accumulation repulsion could significantly enhance 
metastable cluster lifetimes. To the best of our knowl- 
edge, a systematic study of nucleation and growth kinet- 
ics as a function of repulsion strength and range has never 
been performed (see, however, [57]). Such a study would 
be quite useful in understanding the role of the added 
polymer on the stability of the cluster phase observed by 
Lu et al. 



B. Metastable gas-liquid phase separation 

A generic feature of systems with short-ranged at- 
tractions is the existence of a metastable gas-liquid bin- 
odal buried below the gas-crystal coexistence line. This 
feature gives rise to another possible phase separation 
mechanism. In particular, a two-step nucleation process, 
where the first step of phase separation is the formation 
of a higher density liquid and the subsequent step in- 
volves the transformation of the high density liquid into 
a crystal, arises. This two-step nucleation process can 
proceed by a lower free energy pathway than that of 
classical nucleation, thus accelerating the rate of crystal 
formation. At low densities, just below the metastable 
gas-liquid binodal, interesting nonequilibrium behavior 
has been reported. 

In particular, it is in the vicinity of the metastable gas- 
liquid coexistence line that the bead phase of Sedgwick et 
al. is observed [16]. Thus, in this section, we investigate 
the dynamics of domain growth at two volume fractions 
(point C and C in Fig.[3|) just below the metastable gas- 
liquid binodal. 

When our system is quenched to point C of Fig. [3l 
rapid nucleation of ramified liquid regions, followed by 
solidification occurs. Figure [5] shows this two-step pro- 
cess which is clear upon visual inspection, and may also 
be observed in the time evolution of the energy of the sys- 
tem per particle, as shown in Fig. H For lower density 
quenches (point C), liquid beads nucleate quickly. The 
evolution of the low density droplet phase then occurs via 
cluster diffusion. Generically, the droplets begin to crys- 
talize before they coalesce, as seen in Fig. [6l The crystal 
growth process thus takes place by cluster coalescence, 
and it strikingly similar to the growth of the crystal above 
the metastable gas-liquid binodal. It should be noted 
that the entire process illustrated in Fig. [6] occurs spon- 





FIG. 6: (Color online) Fate of nucleated liquid droplets upon 
a rapid quench to point C of Fig. [S] Upper panel shows early 
stages (t = 2.2 x 10^) and lower panel shows later stages of 
solidification process (t = 3.9 x lO'^). For clarity particles 
that are part of small clusters (n < 6) and monomers are 
represented with smaller spheres. 



taneously, in contrast to the ad hoc seeding procedure 
that we have performed to illustrate gas-crystal phase 
separation at higher temperatures above the gas-liquid 
binodal. Indeed, this is possible due to the rapid forma- 
tion of liquid beads in this region of the phase diagram. 
Clearly, we do not observe that liquid beads are long- 
lived along the low temperature side of the metastable 
gas-crystal binodal. This app ears to contrast with the 
results of Sedgwick et al. [1^]. It is interesting to note, 
however, that recent confocal microscopy studies of Lu et 
al. [61] at = 0.05 do show clear coalescence of clusters. 
There are several reasons why the final stages of nucle- 
ation might occur more slowly in the experimental system 
of Sedgwick et al. than in our simulations. First, our sys- 
tem is monodisperse, which facilitates the crystallization 
process. Crystal beads may coalesce more easily due to 
surface faceting, providing a more regular contact area 
between clusters. Second, the lysozyme units of Sedg- 
wick et al. carry some residual charge which may hinder 
crystal formation. Third, sedimentation may play some 
role, and is clearly not modeled in our system. It should 
also be mentioned that the quench rate may play a signif- 
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(a) 

FIG. 7: (Color online) Comparison of solidified structures of LJ systems with n = 50 at = 

line (a) to T = 0.2 at t = 7.3 x 10^ and (b) to T = 0.10 at t = 7.3 x 10^ which correspond, 

Fig. El (c) Evolution of the energy per particle for configurations quenched at points C, D, and E. Note that the configuration 

at point C fully crystallizes. 



15% quenched under the spinodal 
respectively, to points D and E in 



icant role in the formation of such phases, as Sedgwick et 
al. have mentioned. Although we have not made a sys- 
tematic study of the effects of quench rate, our prelimi- 
nary studies show the same behavior for systems rapidly 
quenched to point C and those that reach metastable 
equilibrium just above point C (homogeneous fluid) and 
are then slowly quenched to point C. Thus, we find no 
clear evidence of quench rate dependence in our system 
in this region of the phase diagram. 



C. Deep quenches: gels 

When systems with short-ranged attractions are 
quenched below the metastable gas-liquid binodal curve, 
phase separation characterized by large-scale fluctuations 
sets in. At short times, the system segregates into liquid- 
rich and liquid-poor regions. Since the gas-liquid criti- 
cal point is buried below the gas-solid coexistence line, 
the fluid phase is metastable with respect to the crystal. 
Thus, an additional process, the nucleation of the stable 
crystal phase from the metastable fluid competes with 
spinodal decomposition, and the time dependence of the 
nonequilibrium process is governed by a subtle interplay 
of factors, such as the relative rates of solid nucleation 
within the liquid to the rate of the initial decomposition 
process [26|. If either crystallization or vitrification oc- 
curs locally inside the dense component before the phase 
separation process is complete, the long wavelength seg- 
regation of phases slows down considerably. On time 
scales relevant for colloidal experiments, porous solids 
that bear the imprint of this arrested phase separation 
appear [10, [lil, Hj, E, SI, \hM- Since many recent ex- 
perimental and computational studies of this route to 
colloidal gelation have been carried out, we confine our- 
selves here to some rather qualitative features of this pro- 
cess in this section. The distinction between the past 









FIG. 8: Radial distribution function for system quenched (a) 
to point D (t ~ in gray and at t ^ 5 x 10^ in black) and (b) 
to point E (t- 1 X 10^). 
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work and the work presented here is that a rather pre- 
cise characterization of such behavior with respect to the 
location of various phase boundaries (here the location 
of the metastable gas-liquid binodal) may be made. 

In Fig. [71(a) and Fig. [71(b) two snapshots of the interme- 
diate evolution of the phase separation process are shown. 
These snapshots correspond to direct quenches to points 
D and E in Fig. [3l respectively. While quenches just 
below the metastable gas-liquid binodal allow for facile 
crystallization, the evolution of the system for deeper 
quenches is anomalously slow. This is clearly illustrated 
by the time dependence of the energy per particle as il- 
lustrated in Fig. [71(c). For the intermediate quench (D) 
there is residual evolution of the energy, while for deeper 
quenches (E) near complete arrest is observed. 

Given the weaker bonding relative to temperature at 
point D, it is expected that the system may explore 
deeper metastable states during its slow evolutions. This 
is clearly seen in Fig. [71(c). Correspondingly, a greater 
degree of crystallinity is expected in the interior of the 
porous solid when compared to deeper quenched struc- 
tures that arrest at earlier times during their structural 
evolution. Indeed, as shown in Fig. [8] more resolved crys- 
tal peaks are observed in the radial distribution function 
g{r) for more shallow quenches. 

A topic of great interest of late is the arrest of phase 
separation by vitrification as opposed to crystallization. 
In the investigation here, the monodisperse nature of 
the sample favors crystallization. By making the sam- 
ple polydisperse, arrested phase separation will occur via 
a local glass transition of the dense phase [llj, [45[ . The 
characterization of the precise location of such a glass 
transition is greatly complicated by a variety of features, 
including the explicitly nonequilibrium nature of the pro- 
cess and the fact that the effective density of the com- 
ponent undergoing the glass transition is the thermo- 
dynamically defined bulk density. The use of GEMC 
enabling essentially exact determination of stable and 
metastable phase behavior, perhaps combined with the 
ideas of Ref. [59] might lead to a more precise character- 



ization of the notion and location of a "glass transition" 
under the metastable gas-liquid binodal residing in the 
expanded parameter space that includes quench depen- 
dent and thermodynamic parameters. 



V. CONCLUSIONS 

In this work we have demonstrated the general util- 
ity and feasibility of GEMC for the study of the phase 
behavior of systems with rapidly varying, short-ranged 
attractions. Thus, rather exhaustive study of both the 
equilibrium and nonequilibrium behavior of such systems 
is possible for generic potentials and with precise refer- 
ence for various phase boundaries. Here, we have mainly 
made comparison to the predictions of the SCOZA ap- 
proach and made a preliminary investigation of dynam- 
ics. In light of recent experiments that show interest- 
ing nonequilibrium phases in systems with predominantly 
short-ranged attractions, more work should be performed 
to understand the role of quench rate, polydispersity, 
gravity, hydrodynamics, and other factors that may infiu- 
ence the dynamical behavior of such systems. Another 
avenue worthy of study are systems that possess long- 
ranged repulsion in addition to short-ranged attraction. 
Such systems have been the focus of intense recent study. 
Due to the competing lengthscales in these systems, im- 
plementation of the GEMC method is rather demanding 
for larger systems. A preliminary GEMC investigation of 
gelation and microphase separation in such systems will 
be presented in a future publication [40,] . 
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